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TECHNICAL MEMORANDUM 


ON THE ANALYTICAL DETERMINATION OF RELAXATION MODULUS OF 
VISCOELASTIC MATERIALS BY PRONY'S INTERPOLATION METHOD 

INTRODUCTION 


When analyzing materials with time dependent mechanical properties, it is of 
vital importance that the characterization of such properties be accurately defined. 

For viscoelastic materials, properties such as Relaxation Modulus, Creep Compliance, 
and time -dependent Poisson's Ratio are of utmost importance for numerical and closed- 
form solutions to various problems. 

A popular method of obtaining analytical expressions for these properties con- 
sists of obtaining discrete values as a function of the logarithm of time and curve 
fitting the data to an appropriate expression. Due to the decaying nature of such 
properties as Relaxation Modulus and Poisson’s Ratio, they are conveniently repre- 
sented as a series of exponential functions. A widely used form of these functions is 
the so-called Prony series (due to Gaspard Francois Clair Marie Riche de Prony, 
1755-1839), which can be expressed as, 

n + 

B. e 

i=l 


When using expressions such as this for curve fitting it is noted that there are 
too many unknowns for the amount of equations available for simultaneous solving. 

This generally leads to a trial-and-error approach where values of y. are assumed, 

and an expression for the time- dependent variable is obtained. This process is 
repeated until a satisfactory curve fit has been obtained. 

The intent of this paper is to demonstrate that with the use of the Intergraph 
Interactive Graphics Design System (I.G.D.S.) and its three-dimensional sculptured 
surfaces capability the true Prony method can be efficiently used, and the trial and 
error approach totally eliminated. This is accomplished by curve fitting the test data 
to a third-degree, B -Spline interpolation which can be done quickly and with double 
precision accuracy on the I.G.D.S. The resultant curve can then be divided into the 
desired number of equally spaced increments of time for the application of Prony's 
method. It must be stated at this time that equal time increments are necessary in 
order to apply the method properly. It will also be shown that the accuracy of the 
results is not compromised by the introduction of the intermediate step of the B-spline 
interpolation . 

The method is applied in this paper to a solid rocket propellant (TP-H1148) and 
a fluorocarbon elastomer commonly used as O-ring material (V- 747-75). 

The calculations are performed on a Univac 1108 computer and the graphics on 
a VAX 780 computer, both with double precision accuracy. 


DESCRIPTION OF PRONY'S METHOD 


This method is well documented [7,1,6] and will be explained here as a means 
of describing the Fortran program PRONY which will solve, without trial and error, 
for the exponents y^. 

Let Prony's equation be of the form 



m 

A + B i 
i=l 



( 1 ) 


where k = 0,1,2, ...n; i = 1,2,3, ...m; and n > 2m - 1. We first intend to obtain the 

value of Yj and then the values of A and Bp We must remember at this time that- 

the successive values of time (abscissa) must form an arithmetic progression, that is 

t. = t + kw, where w is the difference between any two successive values of time, 
k o 

Equation (1) can now be expressed as, 



m 

A + E B i 

i=l 


Y i (t o +kw) 

e 


( 2 ) 


It is convenient to use the following abbreviations to simplify equation (2), 


V. = 
i 


YjW 

e 


(3) 


C. = 
i 


B i 


Vo 


(4) 


therefore, equation (1) becomes the set of equations 



m 


= A + £ 

i=l 


C V k 
u i v i 


(5) 


Willers [7] explains how the process of obtaining the y. can be simplified by 

forming the difference between two successive ordinates in order to conveniently 
eliminate A from the approximating equations. As an example, take 4 equidistant 
measurements and 3 exponential terms or k = 0, 1,2,3 and i = 1,2,3; this will lead to 
the following set of equations (5), 
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(6a) 


f (t) ~ A + c i + C 2 + C 3 

o 


f (t) x ~ A + C 1 V 1 + C 2 V 2 + C 3 V 3 


(6b) 


f Ct) 2 = A + C 1 V 1 2 + C 2 V 2 2 + °3 V 3 2 


(6c) 


f (,) 3 = A + C 1 V 1 3 + C 2 V 2 3 + ° 3 V 3 3 


(6d) 


the differences can be expressed as, 


S l-0 = C l (V r 1) + C 2< V 2-» + C 3< V 3- 1 > = f ( t)l • f (t) D 


(7a) 


*2-1 = C l (V r 1,V l + C 2< V 2- 1)V 2 + C 3< V 3- 1)V 3 = f <t) 2 ' f ( t)l 


(7b) 


*3-2 = C^Vj-DVj 2 + C 2 (V 2 -1)V 2 2 + C 3 (V 3 -1)V 3 2 = f (t)j - 


(7c) 


*4-3 = C^Vj-DVj 3 * C 2 (V 2 -1)V 2 3 + 0 3 (V 3 -1)V 3 3 = f (t>4 - f (t , 3 . 


(7d) 


If we let the values of Vj be the roots of the equation 


V ra + s x v m_1 + s 2 V m " 2 + ... + V + s m = 0 


( 8 ) 


then we can solve for the values of y. from equations (3) or 


y. = — log V. 
'i w & e i 


(9) 


When trying to interpolate with three or more exponential terms [m > 3 in equa- 
tion (1)] there is the possibility of encountering complex roots V.. This by no means 

precludes successful curve fitting. It does, however, change the final form of equa- 
tion ( 1) thus forcing the analyst to utilize an expression with a combination of expo- 
nential and harmonic terms [1,7]. On the other hand, if the value of V. is a negative 
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real value then equation (9) becomes unsolvable due to the fact that the natural 
logarithm of a negative number is undefined. Fortunately, it has been found that 
impressive accuracy can be obtained, for the materials in question, by using 2 expo- 
nential terms in equation (1), precluding the necessity for value of m greater than 2. 
Our problem now becomes that of obtaining the values of the S. in equation (8). 
Equations (7) can exist simultaneously if we write 1 


A l-0 S 3 + A 2- 1 S 2 + A 3-2 S 1 + A 4-3 0 


( 10 ) 


If we have a total of n+1 data points and m exponential terms in the Prony series, we 
can obtain a total of n-m+1 equations of the form (10) or 


A. n S + A 9 1 S ,+...+ A ( n . = 0 

1~ 0 m 2“ 1 m-1 ( m+1) -m 


A 2- 1 S m + A 3-2 S m-1 + *" + A (m+2)-(m+l) “ 0 


A 3-2 S m + A 4-3 S m-1 + **• + A (m+3)-(m+2) ~ 0 


( 11 ) 


A (n-m+l)-(n-m) ®m + A (n-m+2)-(n-m+l) ^m-1 + ••• + ^n-(n-l) 


Equations (11) can readily be expressed in matrix form as 


1-0 


2-1 


2-1 


3-2 


‘(m+l)-m 


(m+2)-(m+l) 


‘ (n-m+l)-(n-m) a (n-m+2)-(n-m+l) 


... A 


n-(n-l) 


m 


m 


_A =0 ( 12 ) 


or simply 


[6 ] {S} = 0 


(13a) 


where the 6 matrix contains the coefficients and, in this case, k = l,2,3...n. 

We can express (13a) in tensor notation as 


6 


qr 



(13b) 
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where 


q = 1,2,3 ... n-m+1 
r = 1,2,3 ... m+1 , 

in this manner we have <5 ^ = A^q, 6 22 = A 3 _ 2> etc. From equations (13) we can 

form a set of normal equations [1,7] by the method of least squares which written in 
matrix form are 

[a] {S} = {e} . (14a) 

In the tensor notation we have 

a., S £ = e. (14b) 

where 

j, H = 1,2,3, ... m 


and 



n-m+1 

Z 


q=l 


5 qj V 


E i 


n-m+1 


£ % 

q=l 



( 15 ) 


( 16 ) 


where p = m+1. 

Once the oi^ and e. are obtained, equation (14a) can be solved for {S} by 
premultiplying both sides of the equation by the inverse of [a] or 


(S) = [a]' 1 {e} 


( 17 ) 


Equation (8) can now be solved for the V. and from equation (9) the exponents can 
be obtained. 
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COMPUTER IMPLEMENTATION 


Table 1 is the listing of the FORTRAN program PRONY implementing the Prony 
method. The program reads the x. and y^, in this case time and Relaxation Modulus, 

and forms the a matrix. It is then inverted by means of the Gauss- Jordan reduction 
method and the values of S are obtained as in equation ( 17) . 

Once the S. are obtained, then the Yj can be obtained. It is then a simple step 
more to solve for the constants A and Bj by means of the least-squares method. 

TABLE 1. FORTRAN PROGRAM PRONY 

PEOROB IN 197*ME799( 1 ). PRONY 


1 


PARAMETER N-32 , M-2 , NM 1 *N- 1 , MP 1 = M+ 1 . M2 = 2*M 

2 


PARAMETER NM1M-NM1 -M.M2P2-2*M+2 , M4P4=4*M+4 

3 


DIMENSION dC( 12) 

4 


DOUBLE PRECISION R00T(M2P2 ) , AA ( M2P2 ) , BB ( M2P2 ) 

5 


DOUBLE PRECISION CC ( M4P4 ). W , GAM . GAMMA 

6 


DOUBLE PRECISION V(2).S(MP1) 

7 


DOUBLE PRECISION E(N),T(N),DIFF{ NM1 ) , DELTA ( NM 1M, MP 1 ) 

8 


DOUBLE PRECISION ALPHA ( M .MPI).EPS(M) 

9 


DO 10 K* 1 , N 

10 


READ ( 5 , 100)T ( K ) , E ( K ) 

1 1 


WR I TE { 10, 100)T(K) , E ( K ) 

12 

100 

FORMAT ( ) 

13 

10 

CONTINUE 

14 


DO 11 K-1.NM1 

15 


KP 1 =K+ 1 

16 


DIFF(K)=E(KP1 )-E(K) 

17 

1 1 

CONTINUE 

18 


DO 12 1*1, NM 1M 

19 


DO 12 J-1.MP1 

20 


dd=*I+d-1 ' 

21 


DELTA ( I , d)=DIFF(dd) 

22 

12 

CONTINUE 

23 


DO 14 1-1. M 

24 


DO 14 d* 1 , M 

25 


ALPHA( I, J)-0.0 

26 


00 14 K-1.NM1M 

27 


ALPHA! I , d )* ALPHA! I .d)+DELTA(K, I )*D£LTA(K,d) 

28 

14 

CONTINUE 

29 


DO 15 1-1, M 

30 


EPS! I)-0.0 

31 


00 15 d* 1 , NM 1 M 

32 


EPS! I )*EPS! I )+DELTA(d. I )*DELTA(d.MP1 ) 

33 

15 

CONTINUE 

34 


00 16 I-1.M 

35 


ALPHA! I . MP 1 ) * - 1 .0*EPS! I ) 

36 

16 

CONTINUE 

37 


V! 1 )*4 .0 

38 


CALL OGdR! ALPHA. MP1 .M.M.MP1 .$300, dC. V) 

39 

300 

WRITE ( 6 , 301 )dC! 1 ) 

40 

301 

FORMAT ( ' SYSTEM SOLVED UP TO EQUATION NO. ',13//) 

41 


S< 0=1.0 

42 


00 17 1*1. M 

43 


d*M- 1+2 

44 


S! d )* ALPHA! I , MP 1 ) 

45 

17 

CONTINUE 

46 


00 18 1-1, MP1 

47 


WR I TE ( 6 , 400) I , S! I ) 

48 

400 

FORMAT!' S( ' . 12. ' )-' ,£20.6) 

49 

18 

CONTINUE 

50 


TOL-1 .0- 16 

51 


CALL ROOTZ ! S . M , ROOT , AA , BB . CC , TOL ) 

52 


WRITE! 6,61)! ROOT ( I ) , I- 1 .M2 ) 

53 

61 

FORMAT (20X, 'ROOT' ,5X,2E30.8) 

54 


W*T!3)-T<2) 

55 


DO 19 1-1. M 

56 


NN* 2 * I - 1 

57 


GAM-ROOT ! NN ) 

58 


GAMMA- IDLOG! GAM ) )/W 

59 


WRITE 16.700) GAMMA 

60 

700 

FORMAT! 20X . ' EXPONENT ' . E30. 8 ) 

6 t 

19 

CONTINUE 

62 


STOP 

63 


END 


•BRKPT PRINTS 
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NUMERICAL EXAMPLES 


TP-H 1148 SRM Propellant 

Figure 1 shows a plot of the Relaxation Modulus versus log of time as obtained 
from a stress relaxation test [ 2] . As can be seen , the time range extends from 
1.65958E-3 to 1.65958E+1 min. In order to properly apply this method, this data must 
be plotted in real time and then divided into equal time increments, w. This would be 
close to an impossible task unless performed on a graphics computer with the capabili- 
ties of I.G.D.S. 



LOG OF TIME IN MINUTES 

Figure 1. TP-H 1148 Relaxation Modulus. 

As can be seen in Figure 2a, we have plotted seven of the nine data points 
presented in Figure 1. This is done in order to show with clarity the shape of the 
curve, even though all nine points are used when performing the B- Spline interpola- 
tion shown in Figure 2b. As can be seen, the greatest relaxation occurs during the 
first 0.020 min of the test, and therefore we shall concentrate on this portion of the 
curve . 


By "zooming-in” to the portion of interest, we can take increments of time of 
0.001 min. Figure 3 shows the portion between 0.002 and 0.020 min of Figure 2b. 
The I.G.D.S. can accurately locate the intersection between the B- Spline curve and 
the time increments to yield values of Relaxation Modulus, E(t), at each intersection. 
These values can be read directly off the screen and tabulated as in Table 2. This 
is the basic information required in order to run PRONY . 
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Figure 2b. TP-H1148 B-Spline modulus versus actual time. 





20 



Equation (1) shows the general form of the Prony series used in our discussion. 
For this example let us take m = 2 since at m > 3 we have encountered complex values 
of Vj. By inputing the data from Table 2 into our program we obtain the following 

form of equation (8): 

V 2 - 0.822591 V + 0.095196 = 0 (18) 

where , 

S : = -0.822591 
S 2 = 0.095196 

The solution of equation (18) yields the following roots 
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TABLE 2. TP-H 1148 PROPELLANT RELAXATION MODULUS 


TIME 

E (t) 

(MINUTES) 

(KSI) 

.002 

20.3704 

.003 

18.1835 

.004 

16.0257 

.005 

14.4590 

.006 

13.3989 

.007 

12.6508 

.008 

12.1197 

.009 

11.7508 

.010 

11.5090 

.011 

11.3706 

.012 

11.2952 

.013 

11.2244 

.014 

11.1563 

.015 

11.0909 

.016 

11.0282 

.017 

10.9683 

.018 

10.9112 

.019 

10.8569 

.020 

10.8056 

= 0.6832668 


= 0.1392342 



Substitution of equations (19) into equations (9) yield the values of the required 
exponential term constants. 

Yi = n - ru n- log 0.6832668 = -380.8698 
’ 1 0.001 & e 

Y 2 = 07^1 lo S e °* 1393242 = -1971.8239 


We can now rewrite our Prony equation as 

E(t) = A + Bl e" 380 ’ 8698 1 + B 2 e" 1971 - 8239 1 


(19a) 

(19b) 

(20a) 

(20b) 

( 21 ) 
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where t - time in minutes, and A, Bj_, and B2 are constants to be determined by the 
method of least-squares [5]. These constants are, 


A = 10973.3 

(22a) 

B 1 = 23397.3 

(22b) 

B 2 = -79056.7 

(22c) 


It should be noted that equation (21) is valid only for the time range 0.002 to 0.020 
min and that extrapolation beyond this range can lead to erroneous results. Figure 4 
shows a plot of equation (21) as compared to the B-splined experimental data. As 
can be seen extremely good accuracy is obtained with only 2 terms. 


E (t) 
KSI 



TIME IN MINUTES 

Figure 4. TP-H1148 curve fit (Prony's method). 
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V- 747- 75 Fluorocarbon O-Ring Material 

Figure 5 shows a plot of the Relaxation Modulus versus log of time as obtained 
from a stress relaxation test [4]. Once again the time range extends from 1.65958E-3 
to 1.65958E+1 min. The same procedure as for the TP-H1148 solid propellant is 
followed in order to obtain values of E(t) as a function of equal time increments w. 
Figures 6a and 6b are plots of actual data points and B- Spline interpolation, respec- 
tively . 



LOG OF TIME IN MINUTES 

Figure 5. V- 747- 75 Relaxation Modulus. 

For this example we shall take equal increments of time of 0.0025 min and 
analyze the portion of the curve that extends to 0.0800 min. Figure 7 shows the 
equal time increments, and Table 3 shows their respective values of Relaxation Modulus 
as obtained from the I.G.D.S. With this information we are again ready to run PRONY . 

Again we have found that for m > 3 we encounter complex values of V.; there- 
fore, for m = 2 we obtain the following form of equation (8). 

V 2 - 1.14458 V + 0.255224 = 0 (23) 

with 

5 1 = -1.14458 

5 2 = 0.255224 . 
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/- 747—75 MODULUS VS. ACTUAL TIME IN MINUTES 
FIRST SEVEN DATA POINTS) 



TIME (MINUTES) 

Figure 6b. V-747-75 B-spline modulus versus actual time. 




Figure 7. V- 747- 75 Modulus (first 0.080 min). 


TABLE 3. V-747-75 FLUOROCARBON O-RING MATERIAL RELAXATION 

MODULUS AS OBTAINED FROM FIGURE 7, FOR TIME BETWEEN 
0.0025 AND 0.0800 MINUTES 


TIME E <t) 

(MINUTES) (KSI) 


.0025 3977.9 

0050 3317.1 

0075 3032.1 

.0100 2873.1 

.0125 2772.4 

.0150 2687.4 

.0175 2613.8 

.0200 2551.3 

.0225 2499.9 

.0250 2459.7 

.0275 2430.5 

.0300 2411.6 

.0325 2395.4 

.0350 2379.8 

.0375 2364.8 

.0400 2350.4 

.0425 2336.6 

.0450 2323.6 

.0475 2310.7 

.0500 2298.7 

.0525 2287.2 

.0550 2276.4 

.0575 2266.1 

.0600 2256.4 

.0625 2247.3 

.0650 2238.7 

.0675 2230.6 

.0700 2222.8 

.0725 2215.5 

.0750 2208.1 

.0775 2202.1 

.0800 2196.1 
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The roots of equation (23) are 


V 1 = 0.84113 


(24a) 


V 2 = 0.30345 


(24b) 


The required exponential constants are obtained from equation (9) as 


*1 = 4 '°8e V 1 = OMF IO Be °' 84113 = - 69 ' 204 < 25a > 

v 2 = w lo ^e V 2 = 0.0025 Iog e °' 39345 = "477.017 (25b) 


Our Prony equation for the V-747-75 material during the range 0.0025 < t < 0.0800 
is now 


E 


(t) 


A + B 1 


-69. 204t 
e 


+ B. 


-477.017t 

e 


(26) 


The constants A, B]_, and B2 can be determined by the method of least squares [5], 


2239.88 

(27a) 

= 1247.08 

(27b) 

= 2253.44 

(27c) 


Figure 8 shows a plot of equation (26) as compared to the B-Splined experimental 
data. Again a good correlation between the data and the curve fit is seen. 


DISCUSSION OF RESULTS 


Tables 4 and 5 show the comparison between the test data and the results from 
PRONY . We can see in both cases that there is very good accuracy between test 
points and curve fit points as long as we maintain ourselves between the time range 
under study or as determined by Tables 2 and 3. In both cases the value of E(t) 
departs from the actual values as soon as we leave the time range under study. At 
the beginning of the time range this problem can be overcome by extrapolating the 
the test data on to t = 0 [or t = 0.00001 min to avoid numerical problems when plotting 
Log(t)] . At the end of the time range the accuracy can be maintained by including 
as many time points as required for a specific analysis. 
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Figure 8. V- 747- 75 curve fit (Prony's method) 



TABLE 4. TP-H 1148 RESULT COMPARISON CHART 



Relaxation Modulus 

(ksi) (0.002 < t 

< 0.020) 

Time (min) 

Test Data 

B- Spline 

PRONY 

0.00166 

21.114 

21.114 

20.412 

0.00661 

12.913 

12.913 

12.861 

0.01666 

10.992 

10.992 

11.014 

0.04168 

10.144 

10.144 

10.973 

TABLE 5. 

V- 747- 75 RESULT COMPARISON CHART 

Relaxation Modulus (psi) 
(0.0025 < t <: 0.0800) 

Time (min) 

Test Data 

B- Spline 

PRONY 

0.00166 

4247.00 

4247.00 

2636.05 

0.00611 

3113.00 

3113.00 

3179.15 

0.01666 

2639.00 

2639.00 

2634.39 

0.04168 

2341.00 

2341.00 

2309.58 

0.08318 

2189.00 

2189.00 

2243.82 


Another point of interest is that if any of the roots of equation (8) are negative 
we will not be able to solve equation (9) for that particular y. since the natural 

logarithm of a negative number is undefined. A method to overcome this problem has 
not been addressed by any available source on the implementation of Prony's method 
[1,6,7]. It is the author's belief that due to the nature of equation (9) we must 
require that the values of meet the following requirements : 

1. All V. shall be positive real values or complex conjugates. 

2. If the data to be curve fitted is of a decaying nature as a function of time, 
the V. must be less than 1.00 except for complex conjugate values. 

3. If the data has increasing values of the ordinate with increasing time, then 
the V. must be greater than zero, with at least one greater than 1, except for com- 
plex conjugate values. 
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CONCLUSIONS 


A very good tool for the analytical determination of functions which can be 
represented by a series of exponential terms is available. The method, when utilized 
properly, can provide the analyst accurate expressions of the ordinate in question. 
Although the method leaves the unanswered question of what to do when V. is a real 

negative number, the author feels that enough accuracy for materials such as elas- 
tomers and solid propellants can be obtained to perform an accurate stress and strain 
analysis of these materials under load. 
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APPENDIX A 


CURVE FITTING PRONY'S EQUATION BY THE METHOD OF LEAST SQUARES 


After obtaining the values of the y. exponents, Prony's equation (1) can be 

solved for the A and Bj constants from the data points by using the method of least 
squares [5,7], 

The method, due to Legendre, is applied to Prony's equation as follows. Let 
Prony's equation be 


n 

f(t) = A + £ B. 

i=l 



(A- 1) 


If we want to approximate a curve with f(t) as its discrete data points, we can 
express a deviation for the argument value of t as, 


6 m = f (t) m " f(t) m 
m m m 


(A- 2) 


where m = 1,2,3,... ,k and k is the total number of data points. The sum of the 
squares of all the deviations is expressed as, 


D = Z <«t>m - = Z 


6 


'm "'"'nr * m 

m=l m=l 


(A- 3) 


If we substitute equation (A-l) into equation (A-3), we obtain, 


n 


y.t 

1 t r 


D = Z (A + £ B. e 1 m - f(t) m ) J 
m=l i=l 


(A-4) 


The method of least squares requires that equation (A-4) be a minimum. In order for 
this to be true it is necessary that 


3D _ 3D 
3 A 3Bj 


(A-5) 


From equations (A-5) and (A-4) we can write, 


PASS BLANK NOT FBAW 




BlAW 


b'OT FtLMcO 
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( A-6) 


1 3_D 

2 9 A 


S(A + L B i e Titm -f(t) m ) 3 ' <t, ‘" 
m=l i=l 


9 A 


= 0 


n 


Y,-t. 




(A-7) 


m=l i=l 


where j = 1, 2, 3. . .n. 

Equation (A-7) constitutes a set of n equations or one equation for each expo- 
nential term. These latter two equations can be rewritten as 


k 3f(t)„ Yl t„ 3 f(t) 




m ' 'm 


= L «*>, 


3f(t) 


m 


m=l 


i=l m=l 


3 A " Vk 'm 3 A 

m=l 


(A-8) 


k 3f(t)„ ^ Yi t„ 3f(t) 


m=l ^ i=l m=l 


m 


3 B 


- = z «*>, 9f(t)m 


m=l 


'm 3Bj 


(A-9) 


For the purpose of writing a computer program that will solve for the constants A and 
Bj, equations (A-8) and (A-9) may be expressed in matrix form as follows, 


[C] {A} = {F} 


(A-10) 


where the vector {A} holds the sought unknowns A and Bj and the coefficients of [C] 
and {F} are defined as follows, 


c n = £ 


3f(t) 

Ta“ 


m 


= k 


m=l 




Yi- X t m 3f(t) 


m 


3 A 


m=l 


c - f 9f(t)m 

ll 2^ 3B. 

m=l 1 


(A-ll) 


(A- 12) 


(A- 13) 
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e 


(A- 14) 


c ii=x: 

m=l 


Y- n t 9 f(t) 
'i-l m 


m 


3 B . 


j-1 


For equation (A- 14) only we have 


i = j = 2 , 3 , ... k+1 


whereas for equations (A-12) and (A-13), 


i = 1, 2, 3, . . . k 


The coefficients of the vector {F} are defined as, 


F 1 = L «» 


9 f(t) 


m 


m 3 A 


m=l 


? i = L 


9 f(t) 


m 


m=l 


m 3 B. , 

i- 1 


(A- 15a) 


(A- 15b) 


(A-16) 


(A- 17) 


and in this case the i is defined as in equation (A- 15a). 

Once the C and F coefficients are properly identified, both sides of equation 
(A- 10) can be premultiplied by the inverse of [C] to obtain the solution vector {A} 
or, 

{A} = [C]' 1 {F} (A-18) 

Table A-l shows a printout of the program CFIT applying the method described above. 
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TABLE A- 1. FORTRAN PROGRAM CFIT 


PEOROBIN 197*ME799( 1 ) .CFIT 


1 


PARAMETER K=2,KP1=K+1 , KP2 =K+2 


2 



PARAMETER M=33 , VAL= 100 . , MM= 100 

3 



DIMENSION dC(12),E(M),T(M) 

4 



DOUBLE PRECISION B(K),F(KP1 ) , C ( KP 1 , KP2 ) , 1 

5 



DOUBLE PRECISION V(2),M0D(MM) 

6 



DOUBLE PRECISION TT(MM+1 ) , R(K) 

7 

C 


INITIALIZE MATRIX AND VECTOR 

8 



DO 10 N-1.KP1 

9 



F ( N ) =0 . 0 

10 



DO 10 L= 1 , KP2 

11 



C(N , L ) *0 . 0 

12 


10 

CONTINUE 

13 



R( 1 )*1 18. 159 

14 



R( 2 )=66 1 . 228 

15 

C 


R ( 3 ) =4 1 5 . 608 

16 

c 


R ( 4 ) =500. 

17 



DO 11 I = 1 , M 

18 



READ ( 5 , 100 ) T(I),E(I) 

19 


100 

FORMAT ( ) 

20 


1 1 

CONTINUE 

21 



C( 1 . 1 )=1 .0*M 

22 



DO 20 N=2 , KP 1 

23 



00 20 I=1,M 

24 



C( 1 ,N)=C( 1 , N ) +EXP ( -R(N- 1 ) *T ( I ) ) 

25 



C(N. 1)=C(N, 1 ) + EXP ( -R(N- 1 )*T(I ) ) 

26 


20 

CONTINUE 

27 



DO 30 d = 2 , KP 1 

28 



DO 30 N"2 , KP 1 

29 



DO 30 I=1,M 

30 



C(d,N)“C(J,N)+EXP(-R(N-1 )*T( I ) )*EXP( -R(d 

31 


30 

CONTINUE 

32 



DO 40 I = 1 , M 

33 



F ( 1 ) =F ( 1 ) + E ( I ) 

34 


40 

CONTINUE 

35 



DO 50 N-2 , KP 1 

36 



DO 50 I-1.M 

37 



F(N)=F(N)+E( I ) *EXP( -R(N- 1 )*T( I ) ) 

38 


50 

CONTINUE 

39 



DO 60 N= 1 , KP 1 

40 



C( N , KP2 ) * F (N ) 

41 


60 

CONTINUE 

42 



V( 1 )=4.0 

43 



CALL DGJR(C . KP2 , KP 1 . KP 1 , KP2 , $300 . dC . V ) 

44 


300 

WR I TE ( 6 . 301 )dC( 1 ) 

45 


301 

FORMAT ( ' SYSTEM SOLVED UP TO EQUATION NO 

46 



A~C ( 1 , KP2 ) 

47 



DO 70 N = 2 , KP t 

48 



B(N - 1 )=C(N,KP2 ) 

49 


70 

CONTINUE 

50 



WR I TE ( 6 , 102 )A,B( 1),B(2),B(3) 

51 


102 

FORMAT ( 10X.4E 15.6) 

52 



TMAX=T ( M ) 

53 



TMIN=T ( 1 ) 

54 



DO 80 1=1, MM 

55 



PRONY (I ) *0 . 0 

56 


80 

CONTINUE 

57 



XINC= ( TMAX-TMIN )/VAL 

58 



TT ( 1 )=TMIN 

59 



DO 90 1=1. MM 

60 



DO 101 N-2.KP1 

61 



PRONY ( I ) *PRONY ( I )+B(N- 1 ) *EXP( -R(N- 1 )*TT( 

62 


101 

CONTINUE 

63 



MOD(I)*A+PRONY(I) 

64 



TT ( 1 + 1 )*TT( I ) +XINC 

65 


90 

CONTINUE 

66 



ARG*TMIN 

67 



Y 1 *0.0 

68 



Y2*0. 0 

69 



ARG 1 *0.0 

70 



DO 110 1=1, MM 

71 



WR I TE ( 9 ) ARG1 , Y 1 , Y2 

72 



ARG1*AL0G10( ARG) 

73 



CALL INTRP5 ( ARG ,T,E,M,2,Y1,N) 

74 



Y2=M0D( I ) 

75 



ARG-ARG+XINC 

76 

c 


WR I TE ( 6 , 105 ) ARG , Y 1 , Y2 

77 

c 

105 FORMAT ( 10X , 3E20. 6 ) 

78 


1 10 

CONTINUE 

79 



STOP 

80 



END 


.13//) 
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